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Abstract 

We present numerical results based on a simplified ecological system in evolu- 
tion, showing features of extinction similar to that claimed for the biosystem 
on Earth. In the model each species consists of a population in interaction 
with the others, that reproduces and evolves in time. Each species is simul- 
taneously a predator and a prey in a food chain. Mutations that change the 
interactions are supposed to occur randomly at a low rate. Extinctions of pop- 
ulations result naturally from the predator-prey dynamics. The model is not 

pinned in a fitness variable, and natural selection arises from the dynamics. 
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I. INTRODUCTION 



The evolution of living organisms is a fascinating phenomenon that has intrigued the 
imagination of the scientific and non-scientific community. However, the formulation of 
mathematical models falls necessarily to drastic simplifications. For example, evolution 
has often been considered as a "walk" in a rugged landscape. Following this line, Bak 
and Sneppen (BS) have proposed a model of biological evolution ]]J that has become quite 
interesting to the physics community due to its simplicity and the new insight it provides 
to the problem. It has been shown that this model evolves to a self-organized critical 
state (SOC), and is kept there by the means of avalanches of evolutionary activity. This is 
appealing for a model of biological evolution, since it has been observed that life on Earth 
could be in a SOC state Nevertheless, models based in fitness landscapes, or in a 

concept of fitness different from the biological one, have been criticized from a biological 
point of view ||£|. 

Since one of the characterizing aspects of life, and perhaps the most fundamental one, 
is that of self-replication, it is our belief that more realistic models should involve a dy- 
namic population for each species. The starting point of combining population dynamics 
with evolution is the association of the rates of birth and death and the carrying capacity 
with phenotypes (observable features that arise from the genotype and are, then, subject 
to mutation) [[|. The fitness, namely the expected number of offsprings produced by an 
individual, arise from them. In this way, the process of natural selection is directed by the 
ecological interactions instead of by a non-biological notion of relative fitness. 

Extinction is an essential component of evolution. The great majority of species that 
have ever lived on Earth are now extinct || . There exist competing hypothesis that account 
extinction as originating from within the biosystem, or from external causes -what has been 
called "bad genes or bad luck". In any case, the pattern of extinctions and of surviving 
species or groups of species is certainly an interesting problem to model, to understand, and 
eventually to check with the fossil record. 
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We show in this contribution a simple model of a large ecological system in evolution. 
This produces features of extinction similar to those claimed for the biosystem on Earth. 
We have chosen to study an ecological model in which each species consists of a population 
interacting with the others, that reproduces and evolves in time. The system is supposed 
to be a food chain, and the interactions to be predator-prey. Mutations that change the 
interactions are supposed to occur randomly at a low rate. Extinctions of populations result 
from the predator-prey dynamics. This approach can be thought as middle way between 
the microscopic simulation of "artificial life" by Ray and others [0], and the coarse-grain 
description of models like BS's. 

II. THE MODEL 

Our model ecosystem consists of a number of species that interact and evolve in time. In 
the course of its time evolution the populations grow and shrink following a set of equations. 
Eventually, some of the species become extinct as a result of their interaction with the 
others. Every now and then we change one of the phenotypic features of one of the species, 
mimicking a random mutation of its genome. This modification produces a perturbation in 
the dynamics of the ecosystem, and eventually leads to the extinctions. 

To be more precise, let's consider a simple example of a food web, namely a one- 
dimensional food chain. N species interact in such a way that the species i feeds on the 
species i — 1, and is eaten by the species i + 1. The species 1 is an autotroph: it feeds at 
a constant rate on an "environment". The species N, the top of the chain, is not eaten 
by any species, but dies giving its mass to the environment. Each species has a population 
that evolves in time and interacts with its neighbors in the chain. Furthermore, we consider 
this evolution in discrete time, which is often more realistic than a continuous one || and 
simpler to simulate in a computer. 

As has been said above, each species acts as a predator with respect to the one preceding 
it in the chain, and as a prey with respect to the one following it. As a further simplification, 
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we suppose that there are no intrinsic birth and death rates, apart from those arising from 
the predation and prey contributions. Let's propose the equations governing this behaviors 
||. As a predator, the "population" (a continuous density) of the species i, n\ changes from 
time t to time t + 1 according to: 

An\ = k i n\- 1 n\(l-n\/c l ), (1) 

where ki is a rate of growth of the predator population and q is a carrying capacity that 
accounts for a limitation imposed by the environment. Note that ([!]) includes this carrying 
capacity in a logistic factor to avoid an unbounded growth of the population. Also observe 
that the growth is proportional to the population of preys, without a "satiation" factor. 
Similarly, as a prey, the population of the species i will diminish according to: 

An\ = ^n\n\ +l . (2) 

The parameters fcj, gi and c« are the phenotypic features of our species. In the course of 
the evolution we allow them to change, mimicking random mutations. Moreover, they are 
the same for all the individuals of each single population. We do not model races, traits, 
polymorphisms or any phenotypic variation within a species, and when a mutation occurs 
it is assumed that the whole population "moves" instantly to the new state. In this sense, 
we are modelling the co-evolution of the species and disregarding the evolution of a single 
one as well as other important phenomena like the formation of new species ||. 

Combining the two roles of predator and prey that each species performs, and the special 
status of the ends of the chain, we can write the following set of equations for the evolution 
of the system: 

An\ = kill] n\ (1 — n\/ci) — ginl +1 n l t for i = 1 . . . N 

(3) 

n o = n N+i = 1; 

where we have introduced two fictitious species, and iV + 1, to take account for the border 
condition. 
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We make two simplifications to the system (^): 1) we assume that all the carrying 
capacities are equal, and equal to 1; 2) we assume that ^ = k i+ \. In this way we reduce the 
number of parameters that define the phenotypic features of the ecosystem. 

The dynamics of the system is as follows. At time t — 0, all the populations and 
interactions are chosen at random with uniform distribution in the interval (0,1). Then 
the populations begin to evolve according to the system ([3]). In the course of the evolution 
driven by eq.@ a population can go to zero. As this can happen asymptotically, we consider 
a species extinct if its population drops below a given threshold. This is reasonable since 
actual biological populations are discrete. In order to keep constant the number of species we 
replace an extinct one with a new one, which can be thought as a species coming to occupy 



the niche left by the extinct one |10j. The new population, and the (two) new interactions 
with its neighbors in the chain, are also drawn at random from a uniform distribution in 
the interval (0, 1). On top of this dynamics of predation and extinctions, we introduce 
random mutations. In each time step a mutation is produced with probability p; the species 
to mutate is chosen at random and the mutation itself consists of the replacement of the 
species with a new one, with a new population and new interactions with its predators and 
preys (all random in (0, 1)). 

Observe that we do not introduce the fitness of a species as a dynamical variable. We do 
not even need to compute it from the "phenotypes" fcj. The fitness, the degree of adaptation 
of a species to the ecosystem, arises from the phenotypes, the populations, and the dynamics, 
and it determines whether a species will thrive or become extinct. Chance is introduced by 
the random mutations (and the random replacement of extinct species). It provides the 
material the natural selection works on. This, in turn, determines the survival of the fittest 
by simply eliminating from the system those species that cannot cope with the competing 
environment. We believe in this way we avoid a fundamental problem in the models of 
evolution as a walk in a fitness landscape, namely that the concepts of fitness is not the 
biological one [Q. 
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III. RESULTS OF THE NUMERICAL SIMULATION 



We have run our model for several chain sizes, ranging from 50 to 1000 sites, and for 
times of about 10 7 steps. In the results reported below we let the system evolve, during a 
transient period, from the initial random state to an organized one. 

In fig. [T] we show a typical evolution of the whole population, J2i n% ■ Although each 
population greatly changes in the course of time (what is not shown in the picture), we 
observe that the whole population remains relatively stable. This is due to the saturation 
factor in the predation term of the evolution equations. This whole population shows a 
short time oscillatory dynamics governed by the competition between species through the 
set of equations (§), and a long time evolution characterized by periods of relative stasis 
and periods of fast change. This feature is the effect of mutations and extinction of some 
species. Without the extinctions and mutations, the dynamics of the system should probably 
be chaotic. But it is not this feature that we want to analyze here. Instead, we shall focus 
on the pattern of extinctions. 

As the set of ki represents the phenotypes of the whole ecosystem, its distribution, P(k), 
can be used to characterize its state. Let's observe what happens in the course of time, 
including the transient mentioned above. Initially the ki are chosen at random, and thus its 
distribution is flat in (0, 1), with mean 0.5. This is shown in fig. ^ as a full line. As time 
passes, and as a result of the dynamics, this distribution shifts to a non-uniform one, as 
shown in fig. [| with dashed lines. The whole distribution shifts towards lower values of the 
interaction, showing a tendency of the system to reduce the coupling between the species. 
In the course of the evolution this distribution fluctuates following the pattern of mutations 
and extinctions, but preserves its form. 

Fig. ^| shows the above mentioned fluctuations in P(k) as the evolution of the mean 
value of hi in the system, after the transient. It corresponds to the same run as fig. [I], 
and the same time window is shown. Similarly to that, it displays a pattern of periods of 
stasis interrupted by periods of fast change, but without the short time oscillations displayed 
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by the population. There are periods of stasis of all lengths, to a degree that the unique 
scale of the figure cannot display. This feature of a lack of a typical length will be analyzed 
immediately. Observe in this figure that the mean value oscillates around 0.24, corresponding 
to a distribution like that shown in fig. |2| with a dotted line. 

The extinction events also display this characteristic pattern of periods of stasis and 
periods of change, without a typical size. In order to characterize this, we have chosen the 
time between two consecutive extinction events, r, which distribution is shown in fig. |] for 
several system sizes an probabilities of mutation. Observe that they follow a power law for 
several decades of large values of r, before a region where the effects of the finite size of 
the system start to appear. This is a sign that the system has self-organized into a critical 
state. In other words, the extinction events are distributed in the time axis in such a way 
that the time between extinctions does not have a characteristic duration -as should have if 
the distribution were exponential. Extinctions appear to come in bursts, or avalanches, of 
any size. 

In fig. |5] the pattern of extinction events of the system is seen in the course of time. 
The graphic displays time in the abscissa and the index in the food chain in the ordinate. 
Each dot marks the moment in which a species has become extinct. Each cross, a species 
that suffers a chance mutation. It can be seen that some mutations trigger avalanches of 
extinction, and that these propagate in the "prey" direction. (Bear in mind that an extinct 
species is replaced by a random new one, most probably with a larger population than its 
predecessor, and observe that this has a negative impact in the corresponding prey.) It is 
also apparent that this avalanches have a complex shape in space-time. It is not easy to 
measure their size since, as can be seen in fig. [|, they overlap. See, for example, a mutation 
that is not followed by any avalanche (lower left), another that triggers a very small one 
(lower right), and several that start events of varying size. In any case, let's define a time 
step, At, divide the time axis with it, and count the number of extinctions in each interval. 
Now, let's call the fraction of species that have become extinct in each interval the size, S, 
of the extinction. S will obviously depend on the time step and on the size of the system: 
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S = S(At, N). If the system is in a critical state this function will obey some scaling law on 
the variable N. In fig. |6| we have scaled the distribution of the system size S(N), P(S, N), 
obtained for different system sizes according to the ansatz: 

P(S,N) = N /3 f(S-N v ). (4) 

We can observe that the four curves collapse to a single one, showing the scaling behaviour 
that is typical of a critical state. 

IV. CONCLUSIONS 

We have introduced a simple model of co-evolution and extinction in a food chain. This 
consists of a finite chain of species of predators and preys. Their populations evolve in 
time following Lotka-Volterra-like equations. Evolution is mimicked by randomly changing 
a phenotype. Natural selection is provided by the deterministic behaviour of the dynamical 
system, that produces the extinction of any species that cannot cope with its interactions. 
No relative fitness or fitness landscape had to be invoked. Nevertheless, the pattern of 
extinctions displayed by this toy ecosystem appears to be similar to that proposed for the 
biosystem on Earth. Namely, the system seems to be in a critical state, in which extinctions 
occur in avalanches. The time between extinctions, and the lifetime of any species follow 
distributions that behave like power-laws of time, implying that there is no typical size for 
the time that a species remains in the system. I should be of interest, in a future work, to 
study the precise instability that produces the shift of the distribution of the interactions 
towards low values. The analytical treatment of this problem is currently under study. 

The author greatly acknowledges Ruben Weht for invaluable discussions, and thank Hilda 
Cerdeira for a careful reading of the manuscript. 
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FIGURES 

FIG. 1. Evolution of the whole population of a chain of 100 species. Probability of mutation, 
P-mut = 10 -5 . Observe the superposition of small oscillations and the much larger spikes signaling 
the mutations and extinctions. 

FIG. 2. Distribution, P(k), of the interactions in a chain of 1000 species. Solid line: initial 
distribution; dashed line: t = 5000; dotted line: t = 5x 10 6 . P mu t = 10 -5 . 

FIG. 3. Evolution of the mean value of the interactions, (k), of a chain of 100 species. 
Pmut = 10~ 5 . The plot corresponds to the same run as that of fig. |. 

FIG. 4. The distribution, P(t), of time t between two consecutive extinctions. P m ut = 

FIG. 5. Space-time pattern of extinctions as they occur in a chain of 200 species. The lower 
plot shows a detail of the upper one. Each dot is an extinction event. Crosses, shown by arrows, 
indicate mutations. P mu t = 10 -5 . 

FIG. 6. Scaled distribution of the size of the extinctions, P(S), as a function of the scaled size, 
S. The four curves correspond to systems of N = 50, 100, 200 and 500 species. 
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